*********************************************************************************
* Table 6. Beliefs – Baseline OLS, OLS+IPW for selection, IPWRA (analytic SEs)
*********************************************************************************
clear
set more off
global dir "replication"

        if "$dir"=="ana"  {  
        global root   = "/Users/adiazescobar/Library/CloudStorage/Dropbox/Gender Norms2/Results/"
        global output = "/Users/adiazescobar/Library/CloudStorage/Dropbox/Gender Norms2/Results/"
        }   
        if "$dir"=="ana_imac"  {  
        global root   = "/Users/adaizescobar/Dropbox/Gender Norms2/Results"
        global output = "/Users/adaizescobar/Dropbox/Gender Norms2/Results"
        }   

if "$dir"=="marie"  {
    global root   = "C:\Users\Marie Boltz\Dropbox\Gender Norms2\Results\"
    global output = "C:\Users\Marie Boltz\Dropbox\Gender Norms2\Results\"
}

if "$dir"=="replication" {
    global root   = "."
    global output = "."
    global datadir = "./data"
    global codedir = "./code"
    global outdir  = "./output"
    capture mkdir "$outdir"
}

* ---- QJE-style graph defaults ----
set scheme s2mono
graph set window fontface "Times New Roman"
global GCOMMON graphregion(color(white)) ///
    plotregion(margin(l+8 r+8 t+8 b+8)) ///
    xlabel(, angle(0) labsize(small)) ///
    ylabel(, labsize(small)) ///
    title(, size(small)) subtitle(, size(vsmall)) ///
    note("", size(vsmall))
global GCOMB graphregion(color(white)) imargin(tiny)

****************************************
use "$datadir/dataALLWAVES.dta", clear

run "$codedir/def_samples.do"


global reps 1000

drop if female_w3 != female & female_w3!=.


* Main sample indicator
local sample /*sample3_ext*/ sample4_other

* Belief blocks
local bel1 fob1 sob1 sob1sdb sob1m 
local bel2 fob2 sob2 sob2m sob2p
local bel3 fob3 sob3 sob3sdb sob3m sob3p
local bel4 fob4 fob6
local bel5 fob5 sob5 sob5m sob5p
local bel6 fob7 sob7 sob7m sob7p
local bel7 fob8 sob8 sob8m sob8p
local bel8 e501 e502 e503 e501b e504 e505 e506 e504b
local bel9 fob9

* All belief variables used in selection model
local bels `bel1' `bel2' `bel3' `bel4' `bel6' `bel7'  `bel9'

* Strata covariates (absorbed FE and used in models)
local strat inactiveF employedF fob2M

* Outcomes
local yvars fob2 sob2 sob2m DplurM2 DplurF2 

* Baseline controls (female already here)
local xvars female c1 c2_2_rp c3 c5 f6 age_dif  level_1 level_2 edu_WmoreH inactiveM unemployedM estrato_12 estrato_3 tec_2 

* Diagnostics for selective exposure (who takes D)
reghdfe D `xvars' `bel2' if `sample' == 1, absorb(`strat') vce(cluster id)

* Working control set for main outcome models
local vars_sample3 `xvars' `bel2' 


********************************************************************************
* IPW weights for selection (attrition weights for being in `sample')
********************************************************************************
quietly probit `sample' D `xvars' `bels' `strat', vce(cluster id)
predict ps_attrition, pr 

gen w_`sample' = 1 / ps_attrition
quietly summ `sample'
local p_uncond = r(mean)
gen ws_`sample' = `p_uncond' / ps_attrition   // normalized (if needed)
gen wa_`sample' = 1                           // unweighted

* Default selection weight (you can switch to ws_ if you prefer normalized)
local weight w_`sample'
summ `weight'

local reps $reps

********************************************************************************
* STORAGE MACROS FOR THE THREE TABLES
********************************************************************************
eststo clear
local j = 1

local estlist_ols_base ""   // 6A: Baseline OLS (no weights) – one model per outcome
local mtitles_ols_base ""

local estlist_ols_ipw ""    // 6B: OLS + selection IPW – one model per outcome
local mtitles_ols_ipw ""

local estlist_ipwra ""      // 6C: IPWRA ATT/PO + OLS-IPW D – one model per outcome
local mtitles_ipwra ""

********************************************************************************
* LOOP OVER OUTCOMES
********************************************************************************
foreach y of local yvars {

    /********************************************************************
     * 1. BASELINE OLS (no selection weights): pooled model only
     *    Male/Female enter as stats rows
     ********************************************************************/

    *** (1.1) Pooled OLS – All (unweighted)
    reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1, absorb(`strat') vce(cluster id)
    scalar N_all = e(N)
    eststo OLS_base_`j'

    *** Means of Y by subgroup (unweighted)
    quietly summarize `y'_w3 if `sample' == 1
    scalar ymean_all  = r(mean)

    quietly summarize `y'_w3 if `sample' == 1 & female == 0
    scalar ymean_male = r(mean)

    quietly summarize `y'_w3 if `sample' == 1 & female == 1
    scalar ymean_fem  = r(mean)

    *** (1.2) OLS – Male only (unweighted)
    reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1 & female == 0, ///
        absorb(`strat') vce(cluster id)
    scalar N_male     = e(N)
    scalar b_D_male   = _b[D]
    scalar se_D_male  = _se[D]
    scalar t_D_male   = b_D_male / se_D_male
    scalar p_D_male   = 2*(1 - normal(abs(t_D_male)))

    *** (1.3) OLS – Female only (unweighted)
    reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1 & female == 1, ///
        absorb(`strat') vce(cluster id)
    scalar N_fem      = e(N)
    scalar b_D_fem    = _b[D]
    scalar se_D_fem   = _se[D]
    scalar t_D_fem    = b_D_fem / se_D_fem
    scalar p_D_fem    = 2*(1 - normal(abs(t_D_fem)))

    *** RITEST p-values for D (All, Male, Female), with capture to avoid r(111)
    * All
    capture noisily ritest D _b[D], ///
        cluster(id) strata(`strat') reps(`reps') : ///
        reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1, ///
        absorb(`strat') vce(cluster id)
    scalar p_ri_all = .
    if (_rc == 0) {
        matrix P_all = r(p)
        scalar p_ri_all = round(P_all[1,1], .001)
    }

    * Male
    capture noisily ritest D _b[D], ///
        cluster(id) strata(`strat') reps(`reps') : ///
        reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1 & female == 0, ///
        absorb(`strat') vce(cluster id)
    scalar p_ri_male = .
    if (_rc == 0) {
        matrix P_m = r(p)
        scalar p_ri_male = round(P_m[1,1], .001)
    }

    * Female
    capture noisily ritest D _b[D], ///
        cluster(id) strata(`strat') reps(`reps') : ///
        reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1 & female == 1, ///
        absorb(`strat') vce(cluster id)
    scalar p_ri_fem = .
    if (_rc == 0) {
        matrix P_f = r(p)
        scalar p_ri_fem = round(P_f[1,1], .001)
    }

    *** Attach all scalars to the pooled OLS model for 6A
    estimates restore OLS_base_`j'
    estadd scalar N_all      = N_all,      replace
    estadd scalar N_male     = N_male,     replace
    estadd scalar N_fem      = N_fem,      replace

    estadd scalar ymean_all  = ymean_all,  replace
    estadd scalar ymean_male = ymean_male, replace
    estadd scalar ymean_fem  = ymean_fem,  replace

    estadd scalar b_D_male   = b_D_male,   replace
    estadd scalar se_D_male  = se_D_male,  replace
    estadd scalar p_D_male   = p_D_male,   replace

    estadd scalar b_D_fem    = b_D_fem,    replace
    estadd scalar se_D_fem   = se_D_fem,   replace
    estadd scalar p_D_fem    = p_D_fem,    replace

    estadd scalar p_ri_all   = p_ri_all,   replace
    estadd scalar p_ri_male  = p_ri_male,  replace
    estadd scalar p_ri_fem   = p_ri_fem,   replace

    eststo OLS_base_`j'

    /********************************************************************
     * 2. OLS WITH SELECTION IPW WEIGHTS: pooled model only
     ********************************************************************/

    *** (2.1) Pooled weighted OLS – All
    reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1 [pweight = `weight'], ///
        absorb(`strat') vce(cluster id)
    scalar N_all_w = e(N)
    eststo OLS_ipw_`j'

    *** Weighted means of Y by subgroup
    quietly summarize `y'_w3 [aweight = `weight'] if `sample' == 1
    scalar ymean_all_w  = r(mean)

    quietly summarize `y'_w3 [aweight = `weight'] if `sample' == 1 & female == 0
    scalar ymean_male_w = r(mean)

    quietly summarize `y'_w3 [aweight = `weight'] if `sample' == 1 & female == 1
    scalar ymean_fem_w  = r(mean)

    *** (2.2) Weighted OLS – Male
    reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1 & female == 0 [pweight = `weight'], ///
        absorb(`strat') vce(cluster id)
    scalar N_male_w     = e(N)
    scalar b_D_male_w   = _b[D]
    scalar se_D_male_w  = _se[D]
    scalar t_D_male_w   = b_D_male_w / se_D_male_w
    scalar p_D_male_w   = 2*(1 - normal(abs(t_D_male_w)))

    *** (2.3) Weighted OLS – Female
    reghdfe `y'_w3 D `vars_sample3' `y' ///
        if `sample' == 1 & female == 1 [pweight = `weight'], ///
        absorb(`strat') vce(cluster id)
    scalar N_fem_w      = e(N)
    scalar b_D_fem_w    = _b[D]
    scalar se_D_fem_w   = _se[D]
    scalar t_D_fem_w    = b_D_fem_w / se_D_fem_w
    scalar p_D_fem_w    = 2*(1 - normal(abs(t_D_fem_w)))

    *** Attach scalars to pooled weighted OLS model for 6B
    estimates restore OLS_ipw_`j'
    estadd scalar N_all      = N_all_w,      replace
    estadd scalar N_male     = N_male_w,     replace
    estadd scalar N_fem      = N_fem_w,      replace

    estadd scalar ymean_all  = ymean_all_w,  replace
    estadd scalar ymean_male = ymean_male_w, replace
    estadd scalar ymean_fem  = ymean_fem_w,  replace

    estadd scalar b_D_male   = b_D_male_w,   replace
    estadd scalar se_D_male  = se_D_male_w,  replace
    estadd scalar p_D_male   = p_D_male_w,   replace

    estadd scalar b_D_fem    = b_D_fem_w,    replace
    estadd scalar se_D_fem   = se_D_fem_w,   replace
    estadd scalar p_D_fem    = p_D_fem_w,    replace

    eststo OLS_ipw_`j'

    /********************************************************************
     * 3. IPWRA – ATT/PO for All, Male, Female (analytic SEs)
     *    Use selection weights as pweights; attach to pooled weighted model
     ********************************************************************/

    *** (3.1) IPWRA – All
    teffects ipwra (`y'_w3 `y' `strat') ///
                   (D `vars_sample3' `strat') ///
                   if `sample' == 1 [pweight = `weight'], ///
                   atet vce(cluster id)

    matrix bA = e(b)
    matrix VA = e(V)

    scalar ATET_all    = bA[1,1]
    scalar ATET_se_all = sqrt(VA[1,1])
    scalar ATET_pv_all = 2*(1 - normal(abs(ATET_all/ATET_se_all)))

    scalar PO_all      = bA[1,2]
    scalar PO_se_all   = sqrt(VA[2,2])
    scalar PO_pv_all   = 2*(1 - normal(abs(PO_all/PO_se_all)))

    *** (3.2) IPWRA – Male
    teffects ipwra (`y'_w3 `y' `strat') ///
                   (D `vars_sample3' `strat') ///
                   if `sample' == 1 & female == 0 [pweight = `weight'], ///
                   atet vce(cluster id)

    matrix bA = e(b)
    matrix VA = e(V)

    scalar ATET_male    = bA[1,1]
    scalar ATET_se_male = sqrt(VA[1,1])
    scalar ATET_pv_male = 2*(1 - normal(abs(ATET_male/ATET_se_male)))

    scalar PO_male      = bA[1,2]
    scalar PO_se_male   = sqrt(VA[2,2])
    scalar PO_pv_male   = 2*(1 - normal(abs(PO_male/PO_se_male)))

    *** (3.3) IPWRA – Female
    teffects ipwra (`y'_w3 `y' `strat') ///
                   (D `vars_sample3' `strat') ///
                   if `sample' == 1 & female == 1 [pweight = `weight'], ///
                   atet vce(cluster id)

    matrix bA = e(b)
    matrix VA = e(V)

    scalar ATET_fem    = bA[1,1]
    scalar ATET_se_fem = sqrt(VA[1,1])
    scalar ATET_pv_fem = 2*(1 - normal(abs(ATET_fem/ATET_se_fem)))

    scalar PO_fem      = bA[1,2]
    scalar PO_se_fem   = sqrt(VA[2,2])
    scalar PO_pv_fem   = 2*(1 - normal(abs(PO_fem/PO_se_fem)))

    *** Create IPWRA container: copy pooled weighted OLS and attach ATT/PO scalars
    estimates restore OLS_ipw_`j'
    eststo IPWRA_`j'

    estimates restore IPWRA_`j'
    estadd scalar ATET_all    = ATET_all,    replace
    estadd scalar ATET_se_all = ATET_se_all, replace
    estadd scalar ATET_pv_all = ATET_pv_all, replace

    estadd scalar PO_all      = PO_all,      replace
    estadd scalar PO_se_all   = PO_se_all,   replace
    estadd scalar PO_pv_all   = PO_pv_all,   replace

    estadd scalar ATET_male     = ATET_male,     replace
    estadd scalar ATET_se_male  = ATET_se_male,  replace
    estadd scalar ATET_pv_male  = ATET_pv_male,  replace

    estadd scalar PO_male       = PO_male,       replace
    estadd scalar PO_se_male    = PO_se_male,    replace
    estadd scalar PO_pv_male    = PO_pv_male,    replace

    estadd scalar ATET_fem      = ATET_fem,      replace
    estadd scalar ATET_se_fem   = ATET_se_fem,   replace
    estadd scalar ATET_pv_fem   = ATET_pv_fem,   replace

    estadd scalar PO_fem        = PO_fem,        replace
    estadd scalar PO_se_fem     = PO_se_fem,     replace
    estadd scalar PO_pv_fem     = PO_pv_fem,     replace

    eststo IPWRA_`j'

    ********************************************************************
    * Append models and column titles
    ********************************************************************
    local ylab : variable label `y'_w3
    if "`ylab'" == "" local ylab "`y'_w3"

    local estlist_ols_base "`estlist_ols_base' OLS_base_`j'"
    local mtitles_ols_base `"`mtitles_ols_base' "`ylab'""'

    local estlist_ols_ipw  "`estlist_ols_ipw'  OLS_ipw_`j'"
    local mtitles_ols_ipw  `"`mtitles_ols_ipw'  "`ylab'""'

    local estlist_ipwra    "`estlist_ipwra'    IPWRA_`j'"
    local mtitles_ipwra    `"`mtitles_ipwra'    "`ylab'""'

    local j = `j' + 1
}

********************************************************************************
* EXPORT TABLES
********************************************************************************

*=====================================================
* Table 6A – BASELINE OLS (no selection weights)
*=====================================================
#delimit ;
local opts_ols_base style(tex)  
  varlabels(_cons Constant D "Treatment" female "Female (=1)") 
  cells(b(star fmt(%9.3f)) se(par)) 
  stats(N_all N_male N_fem 
        ymean_all ymean_male ymean_fem
        b_D_male se_D_male p_D_male
        b_D_fem  se_D_fem  p_D_fem
        p_ri_all p_ri_male p_ri_fem,
        fmt(%9.3f)
        label("N (All)" 
              "N (Male)"
              "N (Female)"
              "Mean Dep. Var (All)"
              "Mean Dep. Var (Male)"
              "Mean Dep. Var (Female)"
              "OLS β(D) Male"    "OLS SE Male"    "OLS p-val Male"
              "OLS β(D) Female"  "OLS SE Female"  "OLS p-val Female"
              "RI p-val D (All)" "RI p-val D (Male)" "RI p-val D (Female)")) 
  label collabels(none) msign(--) nolz varwidth(1) modelwidth(14) 
  starlevels(* 0.10 ** 0.05 *** 0.01)
;
#delimit cr

estout `estlist_ols_base' using "$outdir/Table6_Beliefs_OLS_baseline_`sample'.tex", ///
    replace `opts_ols_base' mlabels(`mtitles_ols_base')

*=====================================================
* Table 6B – OLS WITH SELECTION IPW WEIGHTS
*=====================================================
#delimit ;
local opts_ols_ipw style(tex)  
  varlabels(_cons Constant D "Treatment" female "Female (=1)") 
  cells(b(star fmt(%9.3f)) se(par)) 
  stats(N_all N_male N_fem 
        ymean_all ymean_male ymean_fem
        b_D_male se_D_male p_D_male
        b_D_fem  se_D_fem  p_D_fem,
        fmt(%9.3f)
        label("N (All)" 
              "N (Male)"
              "N (Female)"
              "Mean Dep. Var (All, weighted)"
              "Mean Dep. Var (Male, weighted)"
              "Mean Dep. Var (Female, weighted)"
              "OLS-IPW β(D) Male"    "OLS-IPW SE Male"    "OLS-IPW p-val Male"
              "OLS-IPW β(D) Female"  "OLS-IPW SE Female"  "OLS-IPW p-val Female")) 
  label collabels(none) msign(--) nolz varwidth(1) modelwidth(14) 
  starlevels(* 0.10 ** 0.05 *** 0.01)
;
#delimit cr

estout `estlist_ols_ipw' using "$outdir/Table6_Beliefs_OLS_IPW_`sample'.tex", ///
    replace `opts_ols_ipw' mlabels(`mtitles_ols_ipw')

*=====================================================
* Table 6C – IPWRA (analytic SEs), only D in coefficient panel
*=====================================================
#delimit ;
local opts_ipwra style(tex)  
  cells(b(star fmt(%9.3f)) se(par)) 
  stats(ATET_all ATET_se_all ATET_pv_all 
        PO_all   PO_se_all   PO_pv_all
        ATET_male ATET_se_male ATET_pv_male
        PO_male   PO_se_male   PO_pv_male
        ATET_fem  ATET_se_fem  ATET_pv_fem
        PO_fem    PO_se_fem    PO_pv_fem,
        fmt(%9.3f)
        label("ATT All"      "Std. Err. ATT All"   "P-value ATT All"
              "PO All"       "Std. Err. PO All"    "P-value PO All"
              "ATT Male"     "Std. Err. ATT Male"  "P-value ATT Male"
              "PO Male"      "Std. Err. PO Male"   "P-value PO Male"
              "ATT Female"   "Std. Err. ATT Female""P-value ATT Female"
              "PO Female"    "Std. Err. PO Female" "P-value PO Female")) 
  label collabels(none) msign(--) nolz varwidth(1) modelwidth(14) 
  starlevels(* 0.10 ** 0.05 *** 0.01)
;
#delimit cr

estout `estlist_ipwra' using "$outdir/Table6_Beliefs_IPWRA_main_`sample'.tex", ///
    replace `opts_ipwra' mlabels(`mtitles_ipwra') ///
    keep(D) varlabels(D "Treatment") 


